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LONG-TERM  GOALS 

The  overall  technical  goal  is  to  implement  a  1/25°  horizontal  resolution  global  ocean  prediction  system 
based  on  the  HYbrid  Coordinate  Ocean  Model  (HYCOM)  with  tides  and  dynamic  sea  ice.  The 
scientific  goals  include  but  are  not  limited  to  a)  evaluation  of  the  internal  tides  representation  in 
support  of  field  programs,  b)  data  assimilation  in  the  presence  of  tides,  c)  evaluation  of  the  model’s 
ability  to  provide  useful  boundary  conditions  to  high  resolution  coastal  models,  d)  interaction  of  the 
open  ocean  with  ice,  e)  shelf-deep  ocean  interactions,  f)  upper  ocean  physics  including  mixed 
layer/sonic  depth  representation,  g)  mixing,  etc. 

OBJECTIVES 

Work  closely  with  NRL  Stennis  and  NOPP  HYCOM  partners  to  (1)  perform  the  R&D  necessary  to 
develop,  evaluate,  and  investigate  the  dynamics  of  1/25°  global  HYCOM  with  tides  coupled  to  CICE 
(the  Los  Alamos  sea  ice  model)  with  atmospheric  forcing  only,  with  data  assimilation  via  NCODA 
(NRL  Coupled  Ocean  Data  Assimilation),  and  in  forecast  mode  and  (2)  to  incorporate  advances  in 
dynamics  and  physics  from  the  science  community  into  the  HYCOM  code. 

APPROACH 

A  series  of  HYCOM  configurations  is  used  to  (a)  evaluate  internal  tides  representation,  (b)  implement 
a  configuration  for  data  assimilation  in  the  presence  of  tides,  and  (c)  investigate  the  interaction  of  the 
open  ocean  with  ice.  HYCOM  development  is  the  result  of  collaborative  efforts  among  the  Florida 
State  University,  University  of  Miami,  and  the  Naval  Research  Laboratory  (NRL)  as  part  of  the  multi- 
institutional  HYCOM  Consortium  for  Data- Assimilative  Ocean  Modeling  (Bleck,  2002;  Chassignet  et 
al.,  2003;  Halliwell,  2004).  This  effort  was  funded  by  the  National  Ocean  Partnership  Program  (NOPP) 
to  develop  and  evaluate  a  data-assimilative  hybrid  isopycnal-sigma-pressure  (generalized)  coordinate 
ocean  model  (Chassignet  et  ah,  2009).  HYCOM  has  been  configured  globally  and  on  basin  scales  at  up 
to  1/25°  (~3.5  km  mid-latitude)  resolution.  More  details  on  the  latest  global  simulations  can  be  found 
at  http://www.hycom.org  and  in  the  separate  ONR  report  on  NRL  activities  by  A.  Wallcraft. 
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RESULTS 


Internal  wave  representation  in  Oceanic  General  Circulation  Models  (OGCMs) 


In  order  to  investigate  the  impact  of  the  model  grid  spacing  and  vertical  coordinate  choice  on  the 
representation  of  internal  waves  in  OGCMs  and  HYCOM  in  particular,  we  performed  a  multi-model 
study  of  internal  wave  generation,  propagation,  and  evolution.  The  configuration  is  a  stratified  channel 
with  a  barotropic  tide  interacting  with  a  Gaussian  ridge  (Khatiwala,  2000;  Di  Lorenzo  et  al.,  2003; 
Legg  and  Huijst,  2006).  The  models  that  were  evaluated  are  the  HYbrid  Coordinate  Ocean  Model 
(HYCOM,  hybrid  vertical  coordinate),  the  Regional  Ocean  Modeling  System  (ROMS,  terrain¬ 
following  vertical  coordinate),  and  the  MITgcm  (geopotential  vertical  coordinate).  These  models  were 
chosen  because  of  their  different  vertical  coordinates  and  their  wide  usage  by  the  community  for 
global,  regional,  and  coastal  applications.  Figure  1  shows  the  cross-vertical  section  of  the  zonal 
baroclinic  velocity  after  5  days  for  two  of  the  models,  HYCOM  (run  in  isopycnal  mode)  and  ROMS, 
for  different  horizontal  grid  spacing  and  for  a  particular  critical  wave  regime  (the  steepness  parameter, 
ratio  of  the  topography  slope  to  the  internal  wave  beam  angle,  is  close  to  1).  The  model  results  are  in 
reasonably  good  agreement  with  the  theoretical  results  for  high  horizontal  resolution,  but  the  baroclinic 
response  is  considerably  weaker  in  the  coarse  horizontal  resolution  simulations  since  they  are  unable  to 
resolve  the  higher  internal  wave  modes.  In  ROMS  with  the  coarse  horizontal  resolution,  the  grid 
stiffness  ratio  parameter  is  violated  and  smoothing  of  the  bathymetry  is  required.  Such  smoothing  can 
lead  to  an  underestimation  of  the  tidal  conversion  process  up  to  50%  (Di  Lorenzo  et  al.,  2003). 
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Figure  1:  Snapshots  (~4.7  days)  of  cross-vertical  section  of  zonal  baroclinic  velocity  for  HYCOM  (left  panels) 
and  ROMS  (right  panels).  The  horizontal  model  resolution  is  Ax  =1.5  km  for  top  panels,  Ax  =  5  km  for  middle 
panels,  and  Ax  -  10  km  for  bottom  panels.  There  are  100  layers/levels  (Az  =  20  m)  in  the  vertical  for  both 
models. 


Diagnosing  spurious  mixing  in  fixed  coordinate  ocean  models 

Due  to  their  numerical  formulation,  fixed  vertical  coordinate  models  are  not  able  to  preserve  adiabatic 
properties  of  advected  water  parcels  (Griffies  et  al.,  2000).  This  numerically-induced  diapycnal  mixing 
can  in  some  cases  be  of  the  same  magnitude  as  naturally  occurring  mixing.  There  are  only  few  studies 
that  have  attempted  to  document  and  quantify  this  spurious  diapycnal  mixing  (Griffies  et  al.,  2000; 
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Morales  Maqueda,  2007;  Holloway,  2007;  Riemenscheider  and  Legg,  2007;  Burchard  and  Rennau, 
2008).  In  order  to  document  and  quantify  the  spurious  mixing  introduced  by  internal  tides  in  both 
ROMS  and  the  MITgcm,  we  used  two  idealized  configurations:  1)  the  well  known  lock  exchange 
problem  (Haidvogel  and  Beckman,  1999)  as  a  reference  and  2)  the  pure  internal  wave  field  described 
in  the  previous  section.  For  both  experiments,  no  explicit  mixing  is  prescribed  and  thus  the  initial 
water  masses  should  remain  at  all  times.  We  use  the  tracer  flux  method  described  by  Griffies  et  al. 
(2000)  to  compute  the  spurious  numerically-induced  diapycnal  diffusivity  for  both  scenarios.  Figure  2 
and  Figure  3  show  snapshots  of  the  temperature  field  at  t  =  20  hr  for  the  lock  exchange  problem  for 
different  resolution  and  advection  schemes.  The  two  densities  initially  present  have  diffused  and  have 
“leaked”  into  new  density  classes  due  to  the  numerically- induced  mixing.  Figure  4  and  Figure  5  show 
Hovmoller  plots  of  the  diapycnal  diffusivity  using  ROMS  and  the  MITgcm  (respectively)  for  the  same 
experiments.  The  two  figures  show  that  finer  model  grid  resolution  decreases  the  numerical  diffusion, 
but  do  not  particularly  improve  the  under/overshoots  of  temperature  when  high-order  advection 
schemes  (i.e.  U3H  and  3rd  order  upwind  DST)  are  used.  The  periodicity  of  2.5  hr  in  the  diapycnal 
diffusivity  signal  is  due  to  the  natural  seiche  of  the  domain.  These  results  compare  very  well  with  a 
recent  study  by  Burchard  and  Rennau  (2008). 


Figure  6  shows  the  ROMS  tidally  averaged  diapycnal  diffusivity  (tidal  cycle  7)  for  the  internal  wave 
configuration  using  the  MPDATA  scheme.  The  higher  magnitudes  for  diapycnal  diffusivity  are  found 
for  the  densest  density  classes  that  correspond  to  steep  slopes  in  the  ROMS  er-levels.  When  the  model 
resolution  increases,  the  spurious  diapycnal  diffusivities  magnitudes  are  reduced. 
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Figure  2:  ROMS  snapshots  of  the  temperature  field  for  the  3rd  order  upstream  biased  advection  scheme  (left 
panels,  Schepetkin  and  Me  Williams,  2000)  and  the  MPDATA  scheme  (right  panels).  The  top  panels  are  the 
COARSE  experiments  with  Ax  =  2000  m  and  10  cr-levels,  the  middle  panels  are  the  MEDIUM  experiment  with 
Ax  =  500  m  and  40  er-levels,  the  middle  panels  are  the  FINE  experiment  with  Ax  =  125  m  and  160  cr-levels. 
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Figure  3:  MITgcm  snapshots  of  the  temperature  field  for  the  3rd  order  upwind  DST  advection  scheme  (left 
panels)  and  the  Superbee  scheme  (right  panels).  The  top  panels  are  the  COARSE  experiments  withzlx  =  2000  m 
and  10  cr-levels,  the  middle  panels  are  the  MEDIUM  experiment  with  Ax  =  500  m  and  40  cr-levels,  the  middle 
panels  are  the  FINE  experiment  with  Ax  =  125  m  and  160  cr-levels.  The  MITgcm  need  a  higher  viscosity 
parameter  than  ROMS  to  preserve  model  stability. 


Figure  4:  ROMS  Hovmoller  plots  of  the  vertically  averaged  diapycnal  diffusivity  (m2/s).  The  top  panels  are  the 
COARSE  experiments  with  Ax  =  2000  m  and  10  levels,  the  middle  panels  are  the  MEDIUM  experiment  with  Ax 
=  500  m  and  40  levels,  the  middle  panels  are  the  FINE  experiment  with  Ax  =  125  m  and  160  levels. 
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Figure  5:  MITgcm  Hovmoller  plots  of  the  vertically  averaged  diapycnal  diffusivity  (m2/s).  The  top  panels  are 
the  COARSE  experiments  with  Ax  =  2000  m  and  10  levels,  the  middle  panels  are  the  MEDIUM  experiment  with 
Ax  =  500  m  and  40  levels,  the  middle  panels  are  the  FINE  experiment  with  Ax  =  125  m  and  160  levels. 


1024.8 


1025 


1025.2 


1025.4 


1025.6 


Ax  =  1.5  km  Ax  =  5  km 

1024.8 


1025 


1025.2 


1025.4 


0  10  20 
x  10  5(m2/s) 


1025.6 


0  10  20 
x  10~W/s) 


Ax  =  1 0  km 


1024.8 

1025 

1025.2 

1025.4 

1025.6 


0  10  20 
x  10  5(m2/s) 


Figure  6:  Tidally  averaged  (over  tidal  cycle  7)  diapycnal  diffusivity  for  the  ROMS-MPDATA  internal  wave 
experiment  with  100  levels.  Left  panel  is  using  Ax  =  1.5  km,  middle  panel  is  using  Ax  =  5  km,  and  right  panel  is 
using  Ax  =  10  km. 

Data  assimilation  and  tides 

Assimilation  systems  typically  consist  of  multiple  interacting  components  for  data  handling,  pre/post¬ 
processing  module  and  analysis.  The  assimilation  schemes  available  for  HYCOM  were  all  developed 
independently  by  the  different  assimilation  groups  and  each  of  these  schemes  use  different  data 
formats,  file  naming  conventions,  metadata  and  interfaces  between  the  different  components.  Early  on, 
faced  with  the  situation  of  handling  multiple  data  formats  and  interfaces,  we  quickly  realized  that  in 
order  to  be  able  to  compare  these  schemes,  we  needed  to  standardize  their  data  flow  and  data  handling 
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infrastructure.  This  led  to  the  development  of  a  HYCOM  “data  assimilation  toolbox”  (Srinivasan  et  al., 
2009).  This  “toolbox”  is  currently  configured  for  the  Gulf  of  Mexico  (GOM)  which  is  an  ideal  setting 
for  exploring  the  performance  of  assimilation  schemes  because  of  its  circulation  features,  limited  size, 
data  availability  and  interesting  dynamics.  The  GOM  forms  a  semi-enclosed  sea  connected  to  the 
Caribbean  Sea  by  the  Yucatan  Channel  in  the  South,  and  to  the  Atlantic  Ocean  in  the  East  by  the 
Florida  Straits.  It  consists  of  a  deep  central  basin  surrounded  by  a  broad  shelf  and  its  circulation  is 
dominated  by  the  Loop  Current.  The  northern  intrusion  of  the  Loop  Current  into  the  GOM  is  highly 
variable.  At  times,  it  short  circuits  the  GOM  and  turns  east  towards  the  Florida  Straits  soon  after 
entering  the  GOM.  At  other  times,  it  extends  all  the  way  to  the  northeastern  Gulf  of  Mexico 
continental  shelf  forming  a  large  meander.  Eddies  pinch  off  the  meandering  Loop  Current  at  irregular 
intervals,  and  form  anti  cyclonic  warm  core  rings  that  propagate  into  the  western  GOM  where  they 
dissipate.  Due  to  its  limited  size,  the  size  of  the  state  vector  in  high-resolution  GOM-HYCOM 
configurations  is  O  (10  -10  )  and  computationally  expensive  advanced  assimilation  schemes  which 
require  on  the  order  of  O  (50-100)  concurrent  model  runs  can  therefore  be  tested  with  high  horizontal 
resolution  GOM-HYCOM  configurations.  Furthermore,  the  GOM  is  well  sampled  by  both  remotely 
sensed  and  in-situ  data. 

In  order  to  be  able  to  test  the  ability  of  data  assimilation  schemes  to  perform  in  the  presence  of  tides, 
we  added  tides  to  the  “toolbox”.  We  used  a  1/25°  horizontal  configuration  with  20  layers.  The  eastern 
and  southern  boundaries  are  nested  within  the  global  1/12°  HYCOM.  In  order  to  avoid  reflections  due 
to  the  tides  at  the  boundaries,  a  radiative  Flather  condition  is  applied  on  the  barotropic  components  of 
the  velocity.  The  initial  condition  corresponds  to  the  condition  of  August  1st  2008  given  by  one  of  the 
assimilative  simulation  of  the  global  1/12°  configuration  of  HYCOM  (a  strong  loop  current  pattern  is 
present,  Figure  7).  The  atmospheric  forcing  is  given  by  the  Navy  Operational  Global  Atmospheric 
Prediction  System  (NOGAPS).  The  barotropic  tidal  transport  and  elevation  are  prescribed  at  the  open 
boundaries  and  are  derived  from  the  Egbert  and  Erofeeva  (2002)  model.  We  run  the  model  for  30  days 
and  all  the  analysis  are  done  after  day  8,  time  after  which  the  model  has  reached  an  equilibrium  state 
(the  tidally  averaged  energy  in  the  domain  neither  increases  or  decreases).  Figure  7  shows  the  initial 
SSH  for  the  full  considered  domain  and  highlight  the  presence  of  the  loop  current.  We  use  the  t-tide 
software  to  extract  tidal  harmonics.  Figure  8  shows  the  tidal  amplitudes  (colored)  and  tidal  phases 
(contoured)  for:  the  main  semidiurnal  constituent  (M2,  Figure  8a),  the  lunar  diurnal  tidal  constituent 
(Oi,  Figure  8b),  and  a  comparison  to  the  M2  computed  from  GOT99  (resolution  is  1/3°)  (Figure  8c), 
and  the  Oi  computed  from  GOT99  (Figure  8d).  The  location  of  the  M2  amphidromic  point  (zero  tidal 
elevation  at  any  time)  compares  well  with  previous  studies  (Khanta,  2005;  Gouillon  et  al.,  2009)  and 
the  GOT99  (global  ocean  tide  model  from  T/P  altimetry;  Ray,  1999)  data  and  tidal  amplitudes  (~30  cm 
at  the  West  Florida  Shelf).  The  spatial  homogeneity  of  the  Oi  tidal  amplitudes  and  co-tidal  lines  also 
compares  well  to  the  GOT99  data  and  previous  study.  Figure  9  shows  a  snapshot  of  a  cross-vertical 
section  of  meridional  baroclinic  velocity  at  89°W.  A  typical  internal  wave  vertical  structure  can  be 
seen  in  the  center  of  the  basin  where  topography  is  rough  and  where  bottom  velocity  can  reach  up  to 
15  cm/s.  Strong  internal  tides  that  propagates  throughout  the  domain  are  also  generated  at  the  West 
Florida  Shelf  (not  shown  here). 
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Figure  7:  SSH  (m)  of  the  initial  condition  (August  2008)  for  the  considered  domain.  The  gray  solid  line  shows 
the  location  of  the  cross-vertical  section  in  Figure  9. 


Figure  8:  Maps  of  tidal  amplitudes  (colored)  and  phases  (contoured)  for:  a)  the  M2  tidal  constituent,  b)  the  01 
tidal  constituent,  c)  M2  from  GOT99,  and  d)  Ol  from  GOT99.  Amplitudes  are  in  meters  and  phases  in  degrees. 
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Figure  9:  Snapshot  of  a  cross-vertical  section  (89°W)  of  baroclinic  meridional  velocities.  The  abscissa  is  index 
of  the  model  and  ordinate  is  depth  in  meters. 

ARCTIC:  Open  ocean/ice  interaction 

The  major  goal  of  the  work  is  to  develop,  validated,  and  improve  a  fully  coupled  modeling  system  of 
the  Arctic  Ocean  and  sea  ice  and  to  participate  in  the  coordinated  Arctic  Ocean  Model 
Imtercomparison  Project  (AOMIP).  The  ocean  model  component  is  based  on  HYCOM.  The  sea  ice 
component  is  based  on  CICE  4  (Hunke  and  Lipscomb,  2004).  The  model  uses  an  Arctic  dipole  grid  in 
two  configurations:  0.72°  (ARCc  0.72)  and  0.08°(ARCc  0.08)  horizontal  resolution  (courtesy  of 
Wallcraft  and  Posey,  NRL). 

The  initial  configuration  of  the  coarse  resolution  model  (ARCc  0.72)  had  coastline  that  followed  the  50 
m  isobath  with  closed  Bering  Strait  (Figure  10)  with  an  horizontal  resolution  in  the  Arctic  Basin  of 
about  30  km.  The  major  purpose  of  this  configuration  was  to  perform  several  test  runs  to  validate  the 
ocean  -  sea  ice  coupling  and  the  sensitivity  of  the  numerical  solution  to  forcing  at  the  open  boundaries 
(OBs).  A  model  test  simulation  was  run  for  2003-2006  starting  from  PHC  3  climatology 
(, http://psc.apl.washington.edu/POLES ).  At  the  lateral  OBs,  the  model  was  relaxed  to  the 
climatological  fields.  At  the  surface,  the  model  was  forced  with  0.5°  3-hourly  NOGAPS  atmospheric 
fields. 

The  high  resolution  model  (ARCc  0.08)  has  same  domain  configuration  as  ARCc  0.72,  but  the 
coastline  is  at  the  10  m  isobath  and  the  Bering  Strait  is  open.  Horizontal  resolution  in  the  Arctic  is  3  -  5 
km.  A  test  run  was  conducted  for  2005-2008.  The  model  was  initialized  from  1/12°  global  HYCOM 
and  the  lateral  OBs  used  the  fields  from  1/12°  global  HYCOM.  At  the  surface,  the  model  was  forced 
with  0.5°  3-hourly  NOGAPS  atmospheric  fields. 

Both  model  runs  have  been  validated  against  climatological  fields.  We  find  that: 

1)  The  sea  ice  edge  from  simulated  concentration  fields  matches  satellite  derived  ice  edge  fairly  well  in 
the  Arctic  Ocean,  but  that  the  model  tends  to  produce  more  ice  in  the  Greenland  Sea  compared  to 
observations. 

2)  The  ocean  surface  temperature  and  salinity  fields  are  in  agreemet  with  climatological  fields  (PHC  3 
and  GDEM  3)  (Figure  11). 

3)  Both  configurations  fail  to  simulate  the  distribution  and  thermohaline  characteristics  of  the  Atlantic 
water  (defined  as  water  mass  with  T  >  0°C)  in  the  Arctic  Ocean  (Figure  12). 

Possible  reasons  for  misrepresentation  of  the  Atlantic  Layer  in  the  simulations  are  as  follows: 
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(a)  Insufficient  vertical  discretization  in  the  deep  ocean:  Target  densities  are  chosen  such  that  the 
maximum  vertical  resolution  is  given  to  the  ocean  layer  above  the  Atlantic  water. 

(b)  Atlantic  water  inflow  to  the  Arctic  Ocean  is  not  correct. 

(c)  Internal  wave  background  diffusivity  is  to  large. 

(d)  In  the  ARCc  0.72:  unresolved  Siberian  shelf  and  closed  Bering  Strait. 


-5000  -4200  -3400  -2600  -1800  -1000  -200 


Figure  10.  Model  bathymetry  of  the  ARCc  0.72.  Land  is  dark  grey.  Left:  existing  model  domain.  Right:  model 
domain  with  improved  bathymetry:  open  Bering  Strait,  shallow  coastline  (10  m)  resulting  in  a  better 
representation  of  the  Siberian  shelf,  better  resolved  channels  in  the  Canada  Arctic  Archipelago. 


We  are  currently  reconfiguring  ARCc  0.72  using  an  improved  bathymetry  (Figure  1,  right)  and  further 
anlyzing  ARCc  0.08. 


Bathymetry  ARCc0.72  expt  08.0 


Bathymetry  ARCc0.72  expt  02.0 


Figure  11.  Salinity  fields  at  10  m.  Left  and  Middle:  from  ARCc  0.72  and  ARCc  0.08,  respectively,  in  August 
2006.  Right:  Climatological  salinity  field  from  GDEM  3. 
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Figure  12.  Vertical  temperature  distribution  along  the  Arctic  Ocean  -  Greenland  Sea  cross-section  (along  150 
W  and  0W  meridians).  Left:  from  ARCc  0.72.  Note  the  colder  Atlantic  Layer  in  the  Eurasian  Basin  when 
compared  to  climatology.  Right:  from  ARCc  0.08.  The  model  has  failed  to  simulate  any  Atlantic  water  in  the 
Canadian  Basin. 


IMPACT/APPLICATIONS 


The  1/25°  (3.5  km  mid-latitude)  resolution,  first  used  in  some  FY09  global  HYCOM  simulations,  is 
the  highest  so  far  for  a  global  ocean  model  with  high  vertical  resolution.  A  global  ocean  prediction 
system,  based  on  1/25°  global  HYCOM  with  tides,  is  planned  for  real-time  operation  starting  in  2012. 
At  this  resolution,  a  global  ocean  prediction  system  can  directly  provide  boundary  conditions  to  nested 
relocatable  models  with  ~1  km  resolution  anywhere  in  the  world,  a  goal  for  operational  ocean 
prediction  at  NAVOCEANO.  At  present,  regional  and  coastal  models  often  include  barotropic  tidal 
forcing,  but  internal  tides  are  not  included  in  their  open  boundary  conditions.  By  including  tidal 
forcing  and  assimilation  in  a  fully  3-D  global  ocean  model,  an  internal  tide  capability  will  be  provided 
everywhere  and  allow  nested  models  to  include  internal  tides  at  their  open  boundaries. 

TRANSITIONS 

None. 

RELATED  PROJECTS 


The  computational  effort  is  supported  by  DoD  HPC  Challenge  and  non-challenge  grants  of  computer 

time.  In  FY09,  1/25°  and  1/12°  global  HYCOM  ran  under  a  new  FY09-1 1  DoD  HPC  Challenge  grant. 
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